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1.  Introduction 

In  the  Large  Eddy  Simulation  (LES)  of  turbulence,  we  simulate  the  spa- 
tially filtered  velocity  field.  The  governing  equations  are  obtained  by  ap- 
plying a spatial  filter  to  the  Navier-Stokes  equations.  The  set  of  filtered 
equations  are  solved  numerically  to  obtain  the  solution  in  terms  of  filtered 
quantities.  Though  successful  LES  have  been  performed  for  several  differ- 
ent flows  (Lesieur  and  Metais,  1996;  Meneveau  and  Katz,  2000),  LES  for 
wall  bounded  flows  posses  some  critical  challenges  (Piomelli  et  al.,  1989; 
Balaras  et  al.,  1996).  Some  of  the  difficulties  in  LES  modeling  of  wall 
bounded  flows  are  due  to  the  strong  inhomogeneity  of  the  turbulence,  the 
scaling  of  the  largest  scales  near  the  wall,  and  the  strong  filter  inhomogene- 
ity commonly  employed  to  accommodate  the  boundary  conditions.  Because 
of  these  issues,  many  of  the  reported  LES  of  wall  bounded  flows  have  no 
filtering  at  all  in  the  wall  normal  direction.  The  problem  associated  with 
strong  filter  inhomogeneity  can  be  thought  of  as  arising  from  the  incon- 
sistency of  the  requirement  to  represent  a sharp  boundary  in  a simulation 
which  resolves  only  large  scales. 

In  a different  context,  immersed  boundary  methods  are  promising  tools 
for  treating  complex  geometries  in  fluid  dynamics  simulations.  In  these 
methods  the  irregular  boundary  is  embedded  in  a regular  grid  and  a forc- 
ing term  is  included  to  account  for  the  embedded  boundary.  This  forcing  is 
such  that  a prescribed  velocity  is  achieved  at  the  given  surface  boundary. 
The  forcing  is  nonzero  only  in  the  neighborhood  of  the  boundary  (few  grid 
points  near  the  wall),  but  due  to  incompressibility  of  the  flow,  the  forcing 
affects  the  entire  flow  field.  The  forcing  can  be  considered  to  be  singular 
field  in  the  computational  domain.  The  concept  was  implemented  by  Peskin 


389 


390 


A.  DAS  AND  R.  D.  MOSER 


(1972)  to  simulate  flow  around  heart  valves,  and  Goldstein  et.  al.  (1993) 
used  this  technique  to  simulate  flows  in  complex  geometries  using  a pseudo 
spectral  method.  Due  to  the  global  nature  of  the  expansion  functions  used 
in  the  spectral  methods,  discontinuities  produce  Gibbs  phenomenon  in  the 
flow  field.  According  to  Goldstein  et.  al.  (1993),  this  oscillation  does  not 
corrupt  the  flow  field  as  the  flow  evolves  in  time.  Mohd-Yusof  (1997,1998) 
used  a different  expression  for  the  forcing  to  remove  the  severe  time  step 
restriction  of  the  method  developed  by  Goldstein  et.  al.  (1993).  This  ap- 
proach was  used  by  Verzicco  et.  al.  (1998)  to  perform  an  LES  in  complex 
geometries.  However,  in  immersed  boundary  techniques  it  is  in  general  dif- 
ficult to  preserve  higher  order  accuracy  and  high  resolution  as  is  generally 
required  in  turbulent  simulations.  This  difficulty  mainly  arises  due  to  the 
presence  of  discontinuities  at  the  boundaries. 

These  two  related  observations  lead  us  to  propose  a new  approach  for 
embedding  boundaries  in  the  context  of  LES.  In  this  approach,  we  treat  the 
wall  by  ‘filtering  through  it’.  A homogeneous  or  nearly  homogeneous  filter 
is  applied  to  the  flow  field  including  the  embedded  boundary.  As  a result, 
the  boundary  is  no  longer  a sharp  interface,  but  is  diffused  across  the  filter 
width.  To  recover  a Direct  Numerical  Simulation  (DNS),  the  filter  width 
is  made  small  enough  to  resolve  all  relevant  scales  of  motion,  but  remains 
finite  so  that  the  boundary  is  not  a sharp  interface.  In  LES,  the  scale  on 
which  the  boundary  is  resolved  is  consistent  with  the  resolved  scales  of 
motion. 


2.  Governing  equations  and  numerical  methods 

To  illustrate  the  concept  of  ‘filtering  through  the  wall’,  we  consider  the  heat 
equation  on  y G [—1, 1].  In  the  computational  domain,  we  include  a buffer 
region  outside  the  boundaries,  where  the  velocities  are  zero. 


du  _ d2u 
dt  dy 2 


\v\  < 1 


(1) 


u = 0 \y\  > 1 


(2) 


These  two  equations  can  be  combined 

H{y). 


du 

dt 


by  introducing  a Heaviside  function 


— L < y < L 


where 


\y\  < 1 

otherwise 
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A homogeneous  filter  is  applied  on  the  domain  — L to  L,  including  the 
boundary.  For  any  function  /,  the  filtered  function  / is  obtained  by 


/= jL_LG{x-onm 


where  G is  the  filter  kernel.  The  resulting  filtered  equation  is 


d2u  ..  du 

— 2 + G(y  - l)^-(l)  - G(y  + !)«-(-!) 


—L  <y<L 


The  boundary  term,  b is  defined: 


o = G{y- 1)^(1) -G{y+ 1)^-1) 

which  includes  the  unfiltered  derivatives  at  the  boundary.  These  are  not 
known  from  the  simulated  filtered  field.  We  propose  to  estimate  the  stress 
by  applying  a suitable  constraint  on  the  field  in  the  buffer  domain,  \y\  > 1. 
The  boundary  terms  obtained  by  this  technique  are  similar  to  the  forcing 
term  used  by  (Goldstein  et  ah,  1993)  except  the  forcing  is  distributed  in 
the  domain  according  to  the  filter  kernel. 

Now  consider  the  Navier-Stokes  equations  for  incompressible  flows  and 
apply  a similar  technique.  After  filtering  the  equations  and  taking  account 
of  the  filtered  boundary,  the  filtered  equations  are 


dui  _ dui  dp  1 d d _ 

dt  V'3  dxj  dxi  ~ Re  dxj  dxj  i + i 

where  bi  is  the  boundary  term.  This  boundary  term  is  expressed  as 


bi(x ) = / aij(x.')njG(x  - x')  dx! 
JdR 


where  a is  the  stress  at  the  boundary,  including  pressure  and  viscous  stress, 
dR,  is  the  boundary  of  the  fluid  region  R and  rij  is  the  unit  normal  to  the 
surface. 

In  simulating  wall  bounded  turbulent  flow,  a large  number  of  grid  points 
are  required  to  resolve  the  near  wall  layer,  as  a result,  in  many  practical 
applications  it  would  be  difficult  to  apply  an  LES  that  resolves  the  wall 
layer.  In  many  LES  of  wall  bounded  flows,  approximate  boundary  condi- 
tions are  used  to  model  the  effect  of  the  wall  layer  (Balaras  et  al.,  1996). 


392 


A.  DAS  AND  R.  D.  MOSER 


The  approximate  boundary  conditions  are  prescribed  in  terms  of  the  wall 
shear  stress,  so  wall  stresses  must  be  determined  in  terms  of  the  resolved 
velocities. 

In  the  present  formulation,  the  unfiltered  wall  stresses  are  also  required, 
and  for  analogues  reason.  Even  though  the  boundary,  or  the  flow  near  it  is 
not  resolved  in  LES,  it  is  necessary  to  represent  the  transfer  of  momentum 
from  the  flow  to  the  wall,  which  is  governed  by  the  wall  stresses.  In  the 
current  description,  in  which  the  unfiltered  velocity  is  zero  in  the  buffer 
domain,  the  wall  stress  is  the  surface  forcing  required  to  ensure  that  mo- 
mentum is  not  transfered  to  the  buffer  domain.  That  is,  that  the  velocity 
remains  zero.  This  suggests  a technique  for  determining  the  wall  stress. 
Instead  of  defining  a force  to  make  the  velocity  zero  at  the  boundary  as  in 
Mohd-Yusof  (1997),  we  choose  awau  to  minimize  the  transport  of  momen- 
tum to  the  exterior  domain.  To  this  end,  the  wall  stresses  at  each  time  step 
is  defined  by  minimizing 


(4) 


where  the  integral  is  over  the  buffer  domain.  The  |u|2  term  forces  the  energy 

I Ml2 


in  the  buffer  domain  to  be  small,  and  the  a | | term  ensures  that  the 

transfer  of  energy  into  the  domain  is  small.  The  constant  a controls  the 
balance  between  these  two  competing  requirements  and  is  set  to  a value 
of  order  At2.  In  the  Fourier  spectral  method  employed  in  section  3 for  the 
channel  flow,  this  minimization  is  straight  forward  since  it  can  be  done 
independently  for  each  (kx,  kz)  wavenumber,  resulting  in  a 6-parameter 
optimization  in  {oXy,ayy,aZy). 


3.  Results  and  Discussion 

In  this  section,  we  present  the  results  of  two  numerical  tests  of  the  filtered 
boundary  technique. 

In  both  cases,  a plane  channel  flow  is  simulated.  Periodic  boundary 
conditions  are  applied  in  streamwise  (a;)  and  spanwise  (z)  directions,  and 
in  the  wall  normal  direction  ( y ),  a buffer  region  outside  the  wall  is  included, 
and  periodicity  is  imposed  in  the  buffer  region  as  well.  A Fourier  cutoff  filter 
is  applied  in  the  wall  normal  direction  effectively  filtering  the  wall,  and  a 
Fourier  spectral  method  is  used.  A low  storage  second-order  Runga-Kutta 
method  is  used  to  time  discritize  the  nonlinear  terms  and  the  viscous  terms 
are  treated  using  an  integrating  factor. 
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3.1.  THE  EVOLUTION  OF  SMALL  AMPLITUDE  DISTURBANCES 

In  this  test,  evolution  of  a linear  stability  mode  is  illustrated.  The  exact 
solution  in  this  case  has  the  form: 

u{x,  y,  t)  = (1  - y2)  + eRe{ipy  (y)el^ax~ljt^ } , (5) 

v(x,y,t)  = ei?e{ia^(y)el(ax'“wf)},  (6) 

where  ip  is  the  Orr-Sommerfcld  eigenfunction  (appropriately  normalized), 
ui  is  the  complex  frequency  (eigenvalue),  a is  the  prescribed  wave  number, 
and  e is  the  perturbation  amplitude.  The  simulation  was  performed  with 
a = 1 at  Reynolds  numbers  Rec  = 10000.  For  this  case,  the  most  unstable 
mode  has  eigenvalue  ui  = 0.23752649  + i0. 00373967.  An  initial  condition 
was  constructed  from  the  eigenfunction  with  the  most  unstable  mode  and 
then  filtered.  Two  cases  were  run,  with  64  and  256  Fourier  modes  in  the  wall 
normal  direction.  A total  of  20  “points”  are  in  the  buffer  domain  outside 
of  the  channel  (\y\  > 1). 

The  exact  unfiltered  pressure  fluctuations  are  formally  zero  in  the  exte- 
rior, resulting  in  a discontinuity  in  pressure,  and  the  resulting  Gibbs  phe- 
nomenon in  the  filtered  pressure  is  shown  in  figure  lc.  The  wall  normal 
pressure  gradient  appears  in  the  u-momentum  equation,  and  this  quantity 
is  dominated  by  the  filtered  delta  function  at  the  boundary  and  the  result- 
ing Gibbs  phenomenon  (figure  Id).  Yet  the  Gibbs  phenomenon  in  velocity 
perturbations  in  figure  la  and  lb  is  imperceptible.  The  reason  is  that  the 
terms  bv  (figure  le)  has  exactly  the  same  structure  as  the  pressure  gradient 
and  cancel  the  Gibbs  phenomenon  (figure  If).  The  role  of  the  boundary 
terms  in  the  momentum  equation  is  thus  clear.  They  regularize  the  stress 
discontinuities  at  the  wall  (both  pressure  and  viscous  stresses). 

In  these  tests,  the  intent  is  to  use  a filter  width  sufficiently  fine  that  the 
only  input  is  the  filtering  of  the  boundary  (i.e.  a “DNS”). 

The  average  growth  rate  {uj1.  imaginary  part  of  the  complex  frequency) 
of  the  disturbance  was  measured  as  a function  of  time.  For  a = 1,  the 
complex  frequency  (w)  is  given  by 

c jj  = -i(l/v)(dv/dt)  (7) 

Due  to  errors  introduced  in  the  initial  conditions,  there  is  a transient  be- 
fore the  asymptotic  growth  rate  is  reached.  The  solution  with  64  grids  in  y 
direction  under  predicts  the  growth  rate  (.00286091)  and  this  is  apparently 
due  to  inadequate  resolution  in  the  domain.  The  growth  rate  of  the  distur- 
bance is  correctly  reproduced  by  the  simulation  with  Ny  = 256  (.00374844 
vs.  the  value  from  linear  theory  .00373967). 

This  is  a severe  test  for  the  filtered  boundary  approach,  since  the  Fourier 
spectral  method  behaves  poorly  in  presence  of  discontinuities  and  because 
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Figure  1.  Effect  of  the  boundary  terms  in  the  test  of  evolution  of  small  disturbances 
in  a channel  flow,  (a)  filtered  u velocity,  (b)  filtered  v velocity,  (c)  filtered  pressure,  (d) 
pressure  gradient,  (e)  boundary  term  for  v equation,  (f)  pressure  gradient  + boundary 
term.  — real  part, imaginary  part. 


the  growth  rate  is  very  sensitive.  But,  the  boundary  terms  precisely  account 
for  the  poor  behavior  of  the  filter  and  provides  a good  representation  of  the 
filtered  eigenfunction. 

3.2.  TURBULENT  FULLY  DEVELOPED  CHANNEL  FLOW 


To  demonstrate  the  applicability  of  this  technique  in  simulating  turbu- 
lent flow,  a fully  developed  channel  flow  is  computed  on  a 64  x 128  x 64 
grid  with  20  point  in  the  buffer  region.  Based  on  the  wall  shear  velocity  uT 
and  channel  half-width  6,  the  Reynolds  number 
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(a)  (b) 


Figure  2.  (a)Mean  velocity,  (h)«rt„s,  normalized  by  wall  shear  stress.  — present  simu- 
lation Re.T  90, Kim  et.  al.  (1987)  ReT  180. 


of  this  flow  is  92.23.  The  corresponding  Reynolds  number  based  on  the 
centerline  velocity  and  channel  half-width  is  1485.  The  streamwise  and 
spanwise  dimensions  of  the  channel  are  47r<5  and  4/37rd  respectively.  The 
flow  field  was  initialized  by  spatially  filtering  a DNS  flow  field  obtained  at 
ReT  = 180  by  Kim  et  al.  (1987). 

The  governing  equations  were  integrated  until  the  mean  flow  reached 
statistical  equilibrium  (approximately  4 S/nT).  Instantaneous  results  are 
shown  at  non  dimensional  time  7.  The  mean  velocity  distribution  normal- 
ized by  the  wall  shear  velocity  is  shown  in  figure  2a,  along  with  the  data 
reported  by  (Kim  et  al.,  1987)  at  ReT  = 180.  Also  shown  in  figure  2b  is 
the  streamwise  turbulent  intensity.  Note  that  the  computed  peak  is  9.84% 
lower  than  the  Rer  = 180  case,  and  this  is  not  a low  Reynolds  number  effect 
(Keefe  et  al.,  1992).  The  cause  of  this  low  peak  in  urms  is  not  clear  at  this 
point,  but  the  near  wall  resolution  may  be  responsible.  In  other  ways,  the 
turbulence  near  the  wall  is  consistent  with  expectations  for  wall  bounded 
flows.  For  example,  familiar  flow  structures  (streaks,  inclined  shear  layers) 
are  present  near  the  wall  and  are  approximately  of  the  expected  scale. 

4.  Conclusion 

A new  approach  for  treating  boundaries  in  the  context  of  LES,  especially 
extended  boundaries,  was  developed,  and  applied  to  the  computation  of 
turbulent  channel  flow.  In  this  representation,  the  issue  of  highly  inhomo- 
geneous filtering  near  a wall  in  LES  is  obviated,  since  homogeneous  filters 
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can  be  used  since  the  boundary  itself  is  filtered.  In  this  formulation,  the 
unfiltered  stresses  at  the  wall  are  required,  and  are  obtained  from  a mini- 
mization of  perturbations  in  the  buffer  domain. 

The  method  was  successfully  tested  in  the  evolution  of  a Tollmien- 
Schlichting  wave  in  a channel  flow.  Also  a low  Reynolds  number  fully  de- 
veloped turbulent  channel  was  simulated.  The  new  approach  is  promising 
for  the  treatment  of  boundaries  in  LES,  and  the  application  of  optimal 
LES  (Langford  and  Moser,  1999)  to  the  filtered  boundary  formulation  is 
planned. 
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